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The properties of the localized states of a two component Bose-Einstein condensate confined in 
a nonlinear periodic potential [nonlinear optical lattice] are investigated. We reveal the existence 
of new types of solitons and study their stability by means of analytical and numerical approaches. 
The symmetry properties of the localized states with respect to the NOL are also investigated. 
We show that nonlinear optical lattices allow the existence of bright soliton modes with equal 
symmetry in both components, bright localized modes of mixed symmetry type, as well as, dark- 
bright bound states and bright modes on periodic backgrounds. In spite of the quasi ID nature of 
the problem, the fundamental symmetric localized modes undergo a delocalizing transition when the 
strength of the nonlinear optical lattice is varied. This transition is associated with the existence 
of an unstable solution, which exhibits a shrinking (decaying) behavior for slightly overcritical 
(undercritical) variations in the number of atoms. 
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I I. INTRODUCTION 

Bose-Einstein condensates (BEG) in optical lattices 
(OL) have recently attracted a great deal of attention 
due to the possibility of investigating, both at the theo- 
retical and at the experimental level, interesting physical 

' phenomena such as Bloch oscillations. Landau Zener tun- 
neling, Mott transitions, etc. [U, 

\ The interplay between the nonlinearity (intrinsic in 
the interatomic interaction) and the periodic structure 
(induced by the OL) leads to the formation of localized 
states trough the mechanism of the modulational insta- 
bility of the Bloch states at the edges of the Brillouin 
zone of the underlying linear periodic system [^J . These 
states, also known as gap-solitons, can exist in presence 

\ of both attractive and repulsive interactions 0, [a, Q , this 
last fact being possible only due to the presence of the 

■ OL. 

The existence of gap solitons in repulsive BEG was 
\ experimentally demonstrated in Q. The phenomena of 
Bloch oscillations, generation of coherent atomic pulses 
(atom laser) [8j, superfluid-Mott transition ,9^, were also 
experimentally observed. The OLs considered in these 
experiments act as external potentials (and therefore lin- 
early) on the condensate, this introducing an intrinsic 
(state independent) periodicity in the system. In the fol- 
lowing we shall refer to this type of lattices as linear OLs 
(LOLs). In higher dimensions LOLs were shown to be 
very effective in stabilizing localized states against col- 
lapse or decay, leading to the formation of stable multi- 
dimensional solitons [1C) |. 

Besides LOLs, it is also possible to consider nonlin- 
ear OLs (NOLs) with symmetry properties which de- 



pend on the wavefunction characterizing the state of the 
system. A NOL can be obtained by inducing a peri- 
odic spatial variation of the two body interatomic inter- 
action strength (atomic scattering length), leading to a 
periodic space modulation of the nonlinear coefficient in 
the Gross-Pitaevskii equation (GPE) governing the mean 
field dynamics of the ground state. This periodic mod- 
ulation can be experimentally achieved eithe r by means 
of the standard Feshbach resonance method 11], taking 
an external magnetic field near the resonance which is 
spatially periodic [H, [13, [13, [l^, [lB| , or by the optically 
induced Feshbach resonance technique. In the last case 
the nonlinear periodic potential can be produced by two 
counter propagating laser beams with parameters near 
the optically induced Feshbach resonance [13, [11] ■ A pe- 
riodic variation of the laser field intensity in space and a 
proper choice of the resonance detuning lead to a spatial 
dependence of the scattering length ^] and hence to a 
spatial dependent nonlinear coefficient in the GPE. 

Different interesting phenomena occurring in BEG in 
presence of a NOL have already been studied, such as the 
transmission of wave packets through nonlinear barriers, 
generation of atomic solitons, and existence of localized 
states Refs. [13, [11 [H HII • Mathematical properties of 
the ground state and the existence of localized states of 
quasi-lD BEG in NOL have also been studied in [2^.[23j. 
All these studies [2l|, [13, [1^ concern mainly with scalar 
(single component) ID BEG in NOL . The possibility 
of stabilizing multi-dimensional scalar solitons by means 
of NOLs is presently under investigation (preliminary 
studies show that NOLs are unable to stabilize 2D soli- 
tons if the average nonlinearity is negative), while multi- 
component BEGs in NOL have not been considered yet 
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neither theoretically nor experimentally. This last prob- 
lem arises when two or more BEC atomic species interact 
in presence of periodic spatial modulations of the scatter- 
ing lengths, which can occur between the species (inter- 
species) and/or within the species (intra-species). The 
interaction between the two BEC components leads to 
an inter-species NOL which can play a stabilizing role for 
localized states. Spatial modulations of the intra-species 
scattering length (giving rise to intra-species NOLs) can 
also lead to the existence of new types of soliton states. 

The aim of the present paper is to study the properties 
of the localized states of two-component BEC mixtures in 
NOLs. The case of a sinusoidal variation in space of the 
intra- and inter-species scattering lengths will be consid- 
ered. In particular we show the existence of new types of 
solitons and study their stability by means of analytical 
and numerical methods. The symmetry properties of the 
localized with respect to the NOL are also investigated. 
We show that the NOL allows the existence of bright 
soliton modes with equal symmetry in both components, 
bright localized modes of mixed symmetry type, as well 
as bright-dark bound states and bright modes on periodic 
backgrounds. We also show that, in spite of the quasi ID 
nature of the problem, the fundamental symmetric local- 
ized modes undergo a delocalizing transition when the 
strength of the non linear optical lattice is varied. This 
transition is associated with the existence of an unstable 
solution which exhibit a shrinking (decaying) behavior 
for slightly overcritical (undercritical) variations in the 
number of atoms. 

The phenomenon of the delocalizing transition was also 
investigated in [26] for the case of multidimensional sin- 
gle component BEC solitons in LOL and in 25] for the 
case of one-dimensional BEC's with combined linear and 
nonlinear OL's. Delocalizing transitions in binary BEC 
mixtures have not been previously investigated. 

For the analysis of strongly localized modes (i.e. local- 
ized in one or few cells of the NOL) we will apply the vari- 
ational approach which was shown to be effective for such 
type of problems, while for delocalizing transitions and 
broad solitons we use a vectorial Gross-Pitaevskii equa- 
tion averaged over rapid variations in space of the non- 
linear potential. Results are then compared with those 
obtained by direct numerical simulations of the coupled 
CPE system. As numerical tools to investigate the above 
problems we use both self-consistent exact diagonaliza- 
tions |22j] and generalized relaxing methods [28|. 

The paper is organized as follows. In Section II we 
describe the physical model for the two component BEC 
under action of a NOL, based on the optical manipula- 
tion of the scattering length by optically induced Fesh- 
bach resonances. The model equations are introduced in 
the mean field approximation in terms of two coupled ID 
Gross-Pitaevskii equations with intra- and inter-species 
interaction terms. The problem of existence of soliton 
solutions (when the inter- and intra-species atomic scat- 
tering lengths are periodically modulated in space), the 
symmetry properties of localized modes and their stabil- 



ity are discussed in Sections III. The delocalizing transi- 
tions of fundamental modes and the existence of unstable 
solutions associated with them are studied in section IV. 
The analytical predictions are confirmed by direct numer- 
ical simulations of the full CP equation (Sections II- IV). 
Finally, in Section V, the main results of the paper are 
summarized. 



II. THE MODEL 

Two component condensate represent the mixture of 
atoms in the different hyperfine states [1^ [13, HH, [33 ]. 
We consider here the dynamics of two-component BEC 
in presence of a nonlinear optical lattice produced either 
by spatially varying magnetic fields near a Feshbach res- 
onance (FR) value or by optically induced Feshbach reso- 
nances 19]. According to this last approach, the scatter- 
ing length Gs can be optically manipulated if the incident 
light is close to the resonance with one of the bound p 
levels of electronically excited molecules. Virtual radia- 
tive transitions of a pair of interacting atoms to this level, 
can change value and/or reverse the sign of the scattering 
length. The periodic variation of the laser field intensity 
in the standing wave, I{x) = Iq cos'^{kx), produces peri- 
odic variation of the atomic scattering length, such that 



1 
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(1) 



where Gso is the scattering length in the absence of light, 
6 is the frequency detuning of the light from the FR and a 
is a constant factor. For weak intensities, when Iq <C \6\, 
we have that Os = + o,si cos^(fca;). Periodic variation 
of the scattering length by a spatially varying external 
magnetic field B{x) near a Feshbach resonance (FR), can 
be described by 



(2) 



Bo - B{x) 



where Ba is the resonant value and A the corresponding 
width. Examples are: a multicomponent BEC of ^^Na 
atoms 32] or the mixture of ^^K - ®^Rb atoms on the 
surface of a chip. The periodic variation of B can be 
controlled by the current in a magnetic wire on the chip 
surfacejsi]. For the mixture '^^K - ^'^Rb it was shown 
recently that the inter-species scattering length ai2 can 
be tuned using the Feshbach resonances by varying the 
external magnetic field in the interval (50 — 800)G |29l] . 

The mean field equations for the ground state wave- 
function of a quasi-lD two-component BEC under the 
action of a NOL are given by the following coupled CP 
equations [13, [II]: 

ift^ = + 5.(^)1^.1'^. + ^12(X)1^3-.1V.,(3) 

where z = 1, 2 refer to the component index, and the 
full-dimensional space and time variables are given by 
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X and t. In the above, 5*12 is the parameter giving the 
strength of the inter-species NOL and Bi{x) is directly 
related to the atomic scattering length of the species i 
{Bi[x) — 2as^ihu!±). In the following we fix the spatial 
dependence of Bi{x) and Si2{x) as 



B, EE ==r,o + r,cos(2fc5) 

5*12 = Si2{x) ^ Gio + GiCos{2kx). 
Introducing the dimensionless variables 



(4) 



The total energy can be obtained from Eqs. ((G]) and ([7]): 

Ui \ 

^, (9) 



where 
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X = xk,t = tujR,ujR = —,Er = — — ,Ui 

a Zm V £yfi 




and 



Pi EE Pi{x) = + 7j cos(2a;) 
(712 = cri2(x) = g.iO + 5i cos(2a;), 



(5) 



where 7^0 = Tio/lriol, 7» = Ti/lr^ol, 5io = Gio/lTjol, 
giQ = Gi/lFjol. Below we will consider the particular 
case gio = 920 = 5o, 5i = .92- 

We can rewrite the above pair of equations as 
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where the normalization of the total wave-function ^' is 
related to the components itj and to the number of atoms 
Ni by the equation 



dx {ul U2) 



Ul 
U2 



N1 + N2. (7) 



We remark that in experiments the magnitude and sign 
of both the inter- and the intra-species scattering lengths 
can be controlled by external magnetic fields [llj or by 
counter propagating laser fields 0, [sl] . In the case of 
immiscibility, when gi2 < \/\giig22\, the repulsive cross- 
interaction between the components affects strongly the 
self-interaction between components. 
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The corresponding Lagrangian is given by 
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dx 



[lio + li cos(2x))|Mi 



- (50 +.gicos(2a;))|uinu2p. (12) 
In our variational approach we consider ut given by 

[x + {il2~i)xoY^ 



Ui 




exp 



2a? 



(^ = 1,2), 



(13) 

where the normalization Ni is related to the number of 
atoms of the species a.^ is the corresponding width, and 
xo is a parameter given the relative initial position of the 
two components. By substituting this ansatz in Eq. (jlip 
and in the averaged Lagrangian L = Cdx, we obtain: 



-HE 



2a? + 



iViiV2 



G 



1 



■(14) 



A. Variational Approach 
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In this section we perform an analytical study in the 
framework of the variational approach (VA) for the case 
of localized (soliton) solutions of the form: Ui{x,t) = 
Wi(x) exp(— i/iit), where /ii arc the chemical potentials. 
From Eq. ^ we have 
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(50 + 9i cos{2x))ul_^Ui. 
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where 
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From the Euler-Lagrange equations dL/dN = 0, 
dL/da ~ and dL/dxo = we obtain the equations 
for the chemical potentials fii and number of atoms Ni: 
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2£ 



\/tt \/tt 
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(18) 
(19) 



with 
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(20) 



(rj.2 _ al+al 



(2a;otan(a;o^l:p^)-a,4) 



(a? + ai) 



2^2 



[a\ + ai)5/2 
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= ko + (1 + a?)7.ie-"'/2 cos(a;o) 



(21) 
(22) 
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^/?>NlN2{^G / dxo) 
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(23) 



where i 7^ j = 1,2. From HH]) and JTH), it also follows 
that 




(24) 



For the particular choice of parameters for the sym- 
metric case, when 710 = 720 = 7o, 7i =72 = 7, we 
have /ii — fj,j — fi, tti — ttj = a, Ni — Nj = N, and 
1*1=1,2(2;, t) = u±{x,t), with 



u±{x,t) 



N 



: exp 



(x±a:o/2)^ 
2a2 



FIG. 1: Considering the symmetric case with 710 = 720 = — 1, 
71 = 72 = —0.5, go = —1, gi = —1.5, with different values 
of xo, we show the VA results for the number of particles A'^, 
chemical potential and energy E, versus the width a. N is 
given in the upper panel, with fia^ (scale in the Ihs) and Aa^E 
(scale in the rhs) given in the lower panel. 



sin(xo) = ^,(.V2)|£ 
7 ^^0 



(25) 



The equations for the chemical potential /i, energy E, 
and number of atoms TV become {xq ^ 0) 
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2 
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1 N 
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1 
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(26) 
(27) 

{Q = Q,,P=P,) (28) 



TV 



27ra 



7o+5o + (7 + 5i)e""'/'l ,(30) 



.(31) 



[7o + .go + (7 + .9i)e-'''/2(l + «')] 



By using (|3ip in pop for the symmetric case with = 
we obtain 



TV 



27ra 



[a" (70 + 7COs(xo)e-'^'/2(^ ^ ^2^^ ^ 



-3 
2^^2 



(7o+go) + (7 + gi)e-"^/^(l-4) 
(7o+.9o) + (7 + ,9i)e-'^'/2(l + a') 



.(32) 



+ (a' - xl) 



go + gie 2 1 



a2 — 



— ^9. 
e 2^ 



The stability of the soliton solution can be investigated 
by using the Vakhitov-Kolokolov criterion [3^ (in the 
present case, implying that for a stable system we should 
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FIG. 2: VA solution for the chemical potential versus the 
number of particles A^, in the symmetric case with the same 
parameters as in Fig. [1] Exact results are also shown for the 
case of xo = (solid line with empty circles, in both panels). 
Here we observe that the small unstable region, d^/dN > 0, 
presented by the VA is not confirmed by the full numerical 
results. 



FIG. 3: Considering the symmetric case with 710 = 720 = 
— 1, 71 = 72 = —0.5, go = 1, and gi — —1.338926, with 
different values of xo, we show the VA results for the number 
of particles A'^, chemical potential fi and energy E, versus the 
width a. N is given in the upper panel, with /xa^ (scale in the 
Ihs) and 4a^E (scale in the rhs) given in the lower panel. 



have dN/dfi < 0), and also by studying the total energy 
E and chemical potentials fj,i as functions of the width 
a. For the symmetric cases (when Ni = N2 = N and 
fii — [i2 — /i), the results of such study is presented 
in Figs. [1] to m considering an attractive inter-species 
scattering length (go < 0) in Figs. [T] and [2l and repulsive 
inter-species scattering length (30 > 0) in Figs. [3] and HI 
From Figs. [Hand [3] we obtain the behavior of N , chemical 
potential /i and energy E as functions of a. The behavior 
of [i versus iV, in order to check the VK criterion, is shown 
in Figs. OandlH This stability study was done mainly by 
using the variational approach (VA), considering different 
values of the parameter xi^^ which gives the position of 
the soliton solution in respect to the optical lattice. In 
case of xq the VA solutions are also compared with full 
numerical results in Figs. [5] and 2] (solid lines with empty 
circles). As observed, the VA gives a good qualitative 
picture of the exact results, with improved quantitative 
results for large values of |^|. 

The dominant behavior of the chemical potential 
/i and energy E^ as functions of the width a, are removed 
in the bottom panels of Figs. [T] and [3] (by a multiplicative 



factor proportional to a^), in order to enhance their a;o 
dependence. As we can verify, in both the cases, the most 
stable configuration is obtained when = tt. 

As we can see in Fig. [51 the single soliton is stable for 

710 = 720 = -1,71 = 72 = -0.5, go = -1,51 = -1-5. 
The VA predicts the existence of small instability region, 
that is not confirmed by the numerical simulations of the 
system of GP equations. This instability region corre- 
sponds to the broad soliton case with aju > 1, where the 
VA approach is not applicable. 



III. SYMMETRY PROPERTIES OF LOCALIZED 
MODES 

In this section we investigate the symmetry properties 
of localized modes with similar number of particles in 
each component. These modes can be of equal symme- 
try or of mixed symmetry type. In order to find these 
solutions we use both the self-consistent exact diagonal- 
ization method and the generalized relaxing method de- 
scribed in the Appendix (these methods provide identical 
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FIG. 4: VA solution for the chemical potential versus the 
number of particles A^, in the symmetric case with the same 
parameters as in[3] As we can see, the stable region {dfi/dN > 
0) is more pronounced for = tt than for xq = 0. 

results for all the cases studied below, with the only ex- 
ception of the state in Fig. [51 for which the relaxation 
method was not effective). 

In Fig. Owe show the fundamental modes obtained in 
the attractive case (7^0 < 0, 50 < 0) with equal and dif- 
ferent number of atoms in the two components. In both 
cases we have that the maximum of the atomic densities 
are symmetric around the minimum of the corresponding 
effective potentials (see Appendix A). Adopting the same 
terminology introduced in Q for the case of a LOL, 
we shall refer to these modes as OS-OS (onsite symmet- 
ric in both components). Note that OS-OS-modes with 
equal number of atoms have the same chemical poten- 
tials, while for different number of atoms the component 
with a lower number of atoms has also a lower chemical 
potential. For sufficiently strong NOL (see below) these 
modes are very stable under GPE time evolution and rep- 
resent the fundamental ground states of the system in the 
case of all attractive interactions. In particular, the GPE 
time evolution of the density of the OS-OS mode in Fig. [5] 
with a different number of atoms does not show any de- 
viation from the starting density for a time t going from 
to 200. 

Besides onsite symmetry modes, it is also possible to 
have modes that are intersite symmetric (IS) in one or in 



FIG. 5: On-site symmetric modes of Eq. (|6]) with NOL pa- 
rameter values 710 = 720 = —1, 71 = 72 = —0.5, go = — 1, 
gi — —1.5. The top panel refers to the case of an equal 
number of atoms A^i = A^2 = 2 in the two components with 
equal chemical potentials i-ii — I-I2 ~ —3.769. The bot- 
tom panel refers to the case of different number of atoms 
Ni = 2, N2 — 1.5 in the two components and different chem- 
ical potentials /ii — —2.698, fj,2 ~ —2.936. Dashed lines refer 
to second components. 

both components, i.e. symmetric around a maximum of 
the effective nonlinear potential instead than a minimum. 
Such modes can be of type IS-IS (intersite symmetric in 
both components) such as the one shown in the top panel 
of Fig. [S] or of mixed type (OS-IS or IS-OS) such as the 
one shown in the bottom panel of Fig.[6l In contrast with 
the OS-OS mode, the intersite symmetric localized modes 
are found to be unstable under GPE time evolution as 
one can see from Figs. [7] and [5] for IS-IS and OS-IS modes, 
respectively. Notice that in both cases the states decay 
into an OS-OS mode which is the true ground state of 
the system, and that in the IS-OS case the decay of the 
IS component give rise to internal oscillations (relative 
motion between the two final OS components) which can 
last for a long time. Internal oscillations of the OS- OS 
modes can also be excited trough scattering with other 
modes. 

Besides modes that are localized in both components 
it is also possible to couple a localized mode in one com- 
ponent with an extended mode in the other component 
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FIG. 6: Localized modes of the two component GPE in Eq. (|6]) 
of symmetry type IS-IS (top panel) and IS-OS (bottom panel) . 
The number of atoms and chemical potentials of the IS-IS 
mode are A^i = 2, iV2 = 1.5, = -0.653, /i2 = -0.678, 
while for the IS-OS mode are iVi = 0.h,N2 = 1.7, m = 
-0.427, ^^2 = -0.3849. The parameters of the NOL for the 
IS-IS mode are fixed as 710 = 720 = —1, 71 = 72 = 0.5, 
go — —1, gi = 1.5, while for the IS-OS mode are fixed as: 
710 = 720 = -1, 71 = 0.9,72 = -0.5, go = -1.5, gi = 1.5. 
Dashed lines refer to second components. 
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FIG. 7: Time evolution of the modulo square of the IS-IS (top 
panels) and OS-IS mode (bottom panels) depicted in Fig. |6] 
Left panels refer to first components. 



IV. DELOCALIZING TRANSITION OF 
FUNDAMENTAL OS-OS MODES 



such that the extended state acts as a periodic potential 
for the localized mode and forming a bound state. Such 
an example is presented in Figs. [S] and [5] for the case of 
a binary mixture with an average repulsive interaction 
for the first component (710 > 0, I710I > [71 1 ) and an 
average attractive interaction for the second component 
(720 < 0, 1 720 1 > |72|)- This combination of signs for 
the interactions makes the ground state of the system to 
be extended for the first component and locahzed in the 
second one, leading to the formation of the dark-bright 
bound state depicted in the upper panel of Fig. [Sj An- 
other possible solution for the same combination of pa- 
rameters is also verified, as shown in the lower panel of 
Fig. [51 with the formation of a bright-bright state, hav- 
ing one bright solution on top of the background. For 
the considered parameters, both the solutions presented 
in Fig. [5] are quite stable under the GPE time evolution. 
In Fig. [5] we show the time evolution of the dark-bright 
state shown in the upper panel of Fig. [51 



In this section we investigate the existence of a delo- 
calizing transition for the fundamental OS-OS symmetry 
mode. To this regard we recall that for a single compo- 
nent ID BEG with combined LOL and NOL there ex- 
ists a threshold in the number of atoms below which the 
state becomes delocalized. In the limit of rapidly varying 
NOL's one can show, using the averaging method, that 
the system can be reduced to a nonlinear Schrodinger 
equation with cubic and quintic nonlinearities for which 
the existence of a delocalizing transition is known. For 
binary BEG mixtures the same method leads to a coupled 
system of cubic-quintic NLS equations for which delocal- 
izing transitions are also expected to exist. At the transi- 
tion point the localized state becomes spatially more ex- 
tended and displays properties similar to Townes solitons 
of the ID quintic NLS system or of the 2D NLS equa- 
tion with cubic nonlinearity. For broad soliton states, 
i.e. when the soliton width Ig becomes much larger than 
the periodicity scale Ip, we can consider the expansion 
Ui = Ui + eui^i + e^U2,i + with e Ip/h << 1- At the 



8 



-30 



-20 



-10 

X 



10 20 30 



1 .5 



1 .0 



0.5 



0.0 




\Uj(X,t)\ 



\uJx,t)\ 




100 



t 

50 



1 1 

/ 


1 1 


1 














1 






— 1— 








1 , 



-30 -20 -10 
X 



10 20 30 



-30 -20 -10 
X 



10 20 30 



-30 -20 



-10 



10 



20 



30 



FIG. 8; Dark-bright (upper panel) and bright-bright (lower 
panel) soliton modes of Eq. obtained with the NOL pa- 
rameters 710 = 1,720 = —1, 71 = 72 = —0.5, go = —1, and 
gi — —1.5. In both the cases, we have the same chemical 
potentials fii =0.9476 and ^2 ~ —2.746. The resulting num- 
ber of atoms in the two components are, respectively, A'^i =2 
and A'2 =4.5 for the dark-bright solution (upper panel); and 
A^i =64.246 and N2 =0.444 for the bright-bright solution 
(lower panel). 



leading order we obtain 



1 



Ml,, = -cos(2a;)(7,|C/,|"+5i|i7,f )[/,. 



(33) 



Substituting into Eq. ([5]) and averaging over rapid os- 
cillations we get for the slowly varying functions Ui the 
following coupled system with cubic-quintic interactions 

2 2 ^O'l 4 

iUi^t -\- Ui^xx - "yiolUil Uz--go\Uj\ + — -|^7^| Ui + 

O 

^(27. +3i)|[/,|2|i7,f [/, + ^(27, +gi)|t/,|V, = 0. (34) 

When 7io = —go, we obtain a system of coupled quintic 
NLS equations. For the symmetric case Ui = U2 = U, 
the system reduces to the quintic NLS equation 



tUt + U^.,+x\UfU = 0, 



(35) 



with 



X = ItiCti +51) + ^(71 + 3gi + 272). (36) 



FIG. 9: Time evolution of the dark-bright soliton mode de- 
picted in Fig. \8\ 



(37) 



The Townes soliton solution of Eq. ((35|) is 
,3// 1 



C/(M,x) = e'''*(^)i/4- 



X 



with norm given by 



(38) 



This solution behaves as a separatrix between collaps- 
ing and decaying solutions of the quintic NLSE. Here 
dN/dn = and the VK criterion gives marginal stabil- 
ity. The total Hamiltonian is equal zero on this solu- 
tion H{Ut) = 0. For example, for parameters values 
go = — l,.gi = —1.3389,71 = 72 = —0.5, we obtain the 
critical number Nc = 2.41. A comparison with the nu- 
merical results in Fig. [TO] shows that the averaged NLS 
quintic equation overestimates the critical number Nc by 
about 20 percent (notice that for the same parameters val- 
ues we have A^c = 2 in Fig. [TO]) . From this we conclude 
that the quintic NLS can be used only as a qualitative 
model for the delocalizing transitions of two component 
BEC in NOL. In the following we shall investigate Townes 
solitons and delocalizing transitions by recurring to nu- 
merical methods. Delocalizing transitions in binary BEC 
mixtures with NOL and in coupled NLS equations with 
cubic-quintic nonlinearities have not been previously in- 
vestigated. 

To show the existence of this phenomenon in a binary 
BEC mixture in a NOL we vary in time the parameter gi 
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FIG. 10: Top panel. Unstable OS mode of Townes soliton 
type (continuous curve) of Eq. (|6]) at the critical strength 
of the inter-species NOL g\ = —1.3389 for a delocalizing 
transition to occur. The mode has equal critical number of 
atoms A^i = = Nc = 2 with the same chemical potential 
fii = fi2 ~ —0.03498 and same profiles in both components. 
The dashed line represents the corresponding effective poten- 
tial V^^^ = V2^^ in Eqs. (gO]) and (|4l]). Other parameters are 
fixed as 710 = — 1,720 ~ —1, 71 = 72 = —0.5, go — 1. Bottom 
panel. Time evolution of the Townes soliton mode in the top 
panel as obtained from Eq. ©. 



characterizing the intra-species interaction while keeping 
fixed the inter-species NOLs to which the two compo- 
nent are subjected. Starting from a given value of 51 q, 
for which a stable OS-OS mode exists, we adiabatically 
decrease gi to a value giQ — Agi and then increase it back 
to the original value. 

In absence of delocalizing transitions the state will re- 
store to its original form for any decrement Agi, while 
in presence of a delocalizing transition a threshold value 
for Agi will appear above which the state becomes irre- 
versibly delocalized (it cannot be restored to its original 




FIG. 11: Delocalizing transition of an OS-OS mode of a bi- 
nary BEC mixture with equal number of atoms A^i = 2,N2 = 
2. The top panel shows the time evolution of the mode when 
the inter-species parameter gi is varied according to Eq. p9p 
with ti — 50, t2 — 150, Agi = 0.25. The bottom panel shows 
the same evolution but for the case Agi — 0.34. Other pa- 
rameters are fixed as 710 — —1,720 ~ —1, 71 = 72 = —0.5, 
go — 1, giQ — —1.5. The chemical potentials of the initial 
states at t — are fii — fi2 = —0.5401. 

form). In Fig.[Tl]we show the time evolution of an OS-OS 
symmetric state with an equal number of atoms in the 
two components, during a variation of the inter-species 
interaction in time according of the form 

9i{t) = 5io{l - Agi cos[n ^''f'^^'\ }. (39) 

t2 — ti 

From the top panel of this figure it is clear that for 
a small decrement A^i the state is able to restore the 
initial waveform, while for a larger decrement the state 
becomes fully delocalized. In analogy to what has been 
done for the NLS equation with periodic potential and 
quintic nonlinearity [37| . one can characterize the delo- 
calizing transition in terms of the unstable states which 
separate localized modes from extended ones. For the pa- 
rameter used in Fig. 1111 the critical value in the strength 
of the NOL for the occurrence of a delocalizing transition 
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\u.(x,t)\ 




FIG. 12: Time evolution of the Townes soliton mode in Fig.llOl 
for an undercritical Nund = {■999)^Nc (top panel) and over- 
critical Nov = (l.OOl)^A^c (bottom panel) number of atoms. 
Other parameters are fixed as in Fig. 1101 



is found to be Nc = 2. In Fig. [TO] we show the existence 
of an unstable stationary state found in correspondence 
of this value, which has properties similar to the Townes 
solitons of the quintic ID NLS or of cubic multidimen- 
sional NLS. Note that this stationary state corresponds 
to the unstable branch presented in the bottom panel of 
Fig. m [see the exact results for xq — 0]- 

From Fig. [13 it is indeed clear that for slight undercrit- 
ical variations of the norm (number of atoms) the state 
becomes delocalized, while for slight overcritical varia- 
tions of the norm it shrinks to a fully localized mode, 
resembling the behavior of Townes solitons. Notice that 
due to the equal number of atoms A^i = A^2 the modes 
in the two components have identical chemical potentials 
and identical profiles. 

A delocalizing transition is also observed for OS-OS 
states with different number of atoms in the two compo- 




FIG. 13: Onsite symmetric mode of Eq. ^ with unbalanced 
number of atoms A'^i — 2,N2 — 0.5, and for NOL param- 
eters: 7io = -1,720 = -1, 7i = 72 = -0.5, go = -1, 
gi — —1.5. The continuous (dashed) curve refer to the first 
(second) component. The chemical potentials of the modes 
are fii = —1.013, /12 ~ —1.412. The dashed line refers to the 
second component. 



nents. In this case the system shows a much rich behavior 
due to the possibility to use the inter-species interaction 
to stabilize localized states which in absence of interac- 
tion would be extended over the whole system. An ex- 
ample of such inter-species induced localization is given 
in Fig. [13] for an OS-OS symmetric states of Fig. [TO] with 
an unbalanced number of atoms (a large difference in the 
number of atoms in the two components). In particu- 
lar, in absence of the inter-species interactions, the first 
component has enough atoms to be above the delocal- 
izing threshold, while the second component is taken to 
be below such a threshold, so that the state delocalizes 
in absence of interaction. From Fig. [TO] we see, indeed, 
that the presence of the inter-species interaction prevents 
the second component to delocalize, while in absence of 
the inter-species interaction the first component remains 
localized and the second one delocalizes in a quite short 
time. Due to the many parameters of the problem, a 
full investigation of the delocalizing transitions of the 
fundamental OS-OS mode in binary BEG mixtures with 
NOL requires more extensive numerical investigations. 
We plan to do this in a separated publication. 



V. CONCLUSION 

In this paper we have investigated the localized states 
in two-component BEG with periodic modulation in 
space intra-species and inter-species scattering lengths. 
The stability regions are analyzed using the variational 
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FIG. 14: Time evolution of the OS-OS mode in Fig.[T3]in the 
presence (top panels) and in the absence (bottom panels) of 
the inter-species NOL of strength go = —l,gi = —1.5. 



approach and the Vakhitov-Kolokolov criterion. The 
symmetry properties (with respect to the NOL) of the 
localized modes in each component were considered and 
their stability properties investigated. We showed that 
localized modes of OS-OS type are always stable and rep- 
resents the fundamental ground states of the system in 
the presence of attractive interactions. Intersite symmet- 
ric modes and mixed symmetry modes also exist but they 
appear to be metastable under GPE time evolution, de- 
caying into modes of OS-OS-type. The existence regions 
in the parameter space of strongly localized modes (local- 
ized on few cells of the NOL) of fundamental type were 
predicted by mean of the variational ansatz and their sta- 
bility properties predicted by the Vakhitov-Kolokolov cri- 
terion. Localized modes on tops of periodic backgrounds 



and of bright-dark solitons were also shown to exist in 
the case of binary mixtures with opposite interactions in 
the two components. 

In spite of the quasi ID nature of the problem we 
showed that fundamental solitons undergo a delocalizing 
transition when the strength of the intersite non linear 
optical lattice is varied. This transition was associated 
with the existence of an unstable localized solution which 
extends on many lattice cells of the the NOL and which 
exhibit a shrinking (decaying) behavior for slightly over- 
critical (undercritical) variations in the number of atoms. 

This behavior was shown to exists for fundamental 
modes both with equal and unequal numbers of atoms 
in the two components. 

The existence of the delocalizing transition for the fun- 
damental modes was inferred also from a reduced vec- 
tor GPE obtained by averaging the original GPE sys- 
tem with respect to the rapid spatial oscillations intro- 
duced by the NOL. The process of averaging the NOL in- 
troduces high order nonlinearities (cubic-quintic) which 
make the problem to be effectively equivalent to an higher 
dimensional vector GPE system for which delocalizing 
transition, in analogy to single component multidimen- 
sional cases, are usually expected. 

The study of the delocalizing transition for fundamen- 
tal multi-component solitons in terms of an averaged vec- 
tor GPE with higher order nonlinearities, as well as the 
extension of the above analysis to the multidimensional 
case, appear to be interesting problems which deserve 
further investigations. 
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Appendix - Numerical Approach 

The numerical methods employed in this paper are de- 
scribed in the following subsections. 



A. Self-consistent diagonalization algorithm 

We solve the nonlinear eigenvalue problem in ^ by 
treating the nonlinear part in self consistent manner. 
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This amounts to consider the following linear eigenvalue 
problem 



liiui 



d^U2 



(40) 
(41) 



with the effective potentials defined as V^^^ = 
^i{x)\ui\^ + g{x)\u^-i\^ , i = 1,2. To solve these eigen- 
value problem we adopt a discrete variable representa- 
tion [ill and diagonalize the operators Hi = K + V^^^ in 
the discrete coordinate space representation {x„ = no}, 
n = 1,...,-/Vp, a = L/Np. Here K denotes the kinetic 

energy operator K = -^-§^: L is the length of the system 
and Np the number of grid points. By taking as a basis 
the set of vectors |a;„ >= (0, ...0, 1, 0, ...0), n — 1, Np 
and noting that V^'^'^ is already diagonal in this basis 
while K is diagonal in the momentum representation 

of Hi can written as 



fc„<5n,m, we have that the matrix elements 



(XnmXr,,) = {Xn\F-^kF\Xra) + V^^ ^ {na)5n^m, (42) 

where F\xn) denotes the Fourier (unitary) transform of 
the vector |a;„). Standard diagonalization routines are 
then used to find eigenvalues (chemical potentials) and 
eigenf unctions. The nonlinear eigenvalue is then solved 
in self-consistent manner starting from trial wavefunc- 
tions ui, U2, calculating the effective potentials , solving 
the eigenvalue problems (|40p by diagonalizing the corre- 
spondig matrices ((i^ . selecting given eigenstates as new 
trial functions, and iterating the procedure until conver- 
gence is reached (see Refs.^3)lll] for applications to sin- 
gle and multicomponent BEC cases). 



B. Relaxation technique 

The method of relaxation technique was used to check 
the results obtained with the previous method, to im- 
prove their accuracy, and also to make a complete study 
on the stability of the solutions. 

Stable states are obtained using standard relaxation 
algorithm in imaginary time propagation, fixing the nor- 
malizations given by number of atoms of the two species, 
iVi and N2, and obtaining the chemical potentials /ii and 
/X2. For the hyperbolic (unstable) states we extended 
to a coupled equation system the method developed in 
Ref . [28j . scheme C, in which the idea of "back renormal- 
ization" was used. In this method, it is given the chemical 
potential to obtain the number of atoms. 

For a coupled system, the scheme C of Ref. [11] can 
be generalized, evolving the following equations in imag- 
inary time:' 



d(p 



+ Nil3i\ipf + N2ai2\ 



Ml f (43) 



dr \ dx'^ 

where we have normalized and (j) to one, such that 
tp = ui/^/Ni and cj) = U2/y/N2. f3i and ai2 are given by 
Eq. ©. 

In discretized version, the coupled equations (|43p and 
takes the form 



n+l/3 
n+2/3 

i-l-1/3 
1-1-2/3 
An+2/3 



Ar 



Tl+l/3 



OcNf 

^«+2/3 + ^ _ /3i7V{>"|' - ^12N^W'\^) 



Ar . 

— (M2 



/32A^2"l'^"l'-^i2^i"l¥'"l') r 



o 



CNf 



-1/3 



n+l 



n+l 



,n+l 



dx\(j)' 



n+ll2 



where the superscripts {n, n-\-l, etc) refer to time steps. 
OcN is the Crank-Nicolson evolution operation corre- 
sponding to —d^/dx^. Note that, in this coupled system 
the back renormalization (of iV"^^ and A''^^^) is done by 
exchanging the corresponding wavefunctions (as A'^i is as- 
sociated to Lp and A'^2 to 0). This procedure is required 
for stability, as verified in numerical tests. 

The excited states IS-OS and IS-IS depicted in Fig. [H] 
can be obtained by relaxing Eqs. (I43 |) - (l44|) for a; > and 
imposing the Von Neumann boundary conditions in the 
origin, i.e., at x = 0, dip/dx = and d(j)/dx = 0. The 
present relaxation algorithms are unable to find the state 
shown in Fig. [51 which was obtained by the approach 
given in subsection A. 

As compared to the scheme shown in subsection A, the 
advantage of relaxation methods relies on the possibility 
of generalization to higher dimensions with few compu- 
tational resources. 
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